#include <iostream>
#include <vector>
#include <cmath>
#include <algorithm>
#include <numeric>

std::vector<double> generateAnalogSignal(double freq, double duration, double sampleRate) {
    std::vector<double> signal;
    int totalSamples = static_cast<int>(duration * sampleRate);
    for (int i = 0; i < totalSamples; ++i) {
        double t = i / sampleRate;
        signal.push_back(sin(2 * M_PI * freq * t));
    }
    return signal;
}

std::vector<double> sampleSignal(const std::vector<double>& analogSignal, double analogSampleRate, double digitalSampleRate) {
    std::vector<double> sampledSignal;
    int step = static_cast<int>(analogSampleRate / digitalSampleRate);
    for (size_t i = 0; i < analogSignal.size(); i += step) {
        sampledSignal.push_back(analogSignal[i]);
    }
    return sampledSignal;
}

std::vector<double> reconstructSignal(const std::vector<double>& sampledSignal, double originalSampleRate, double targetSampleRate) {
    std::vector<double> reconstructed;
    int interpolationFactor = static_cast<int>(targetSampleRate / originalSampleRate);
    
    for (size_t i = 0; i < sampledSignal.size(); ++i) {
        for (int j = 0; j < interpolationFactor; ++j) {
            reconstructed.push_back(sampledSignal[i]);
        }
    }
    return reconstructed;
}

double calculateError(const std::vector<double>& original, const std::vector<double>& reconstructed) {
    double error = 0.0;
    size_t minSize = std::min(original.size(), reconstructed.size());
    for (size_t i = 0; i < minSize; ++i) {
        error += pow(original[i] - reconstructed[i], 2);
    }
    return error / minSize;
}

int main() {
    double signalFreq = 5.0;
    double duration = 1.0;
    double analogSampleRate = 1000.0;

    auto originalSignal = generateAnalogSignal(signalFreq, duration, analogSampleRate);

    double nyquistSampleRate = 2.1 * signalFreq;
    auto sampledNyquist = sampleSignal(originalSignal, analogSampleRate, nyquistSampleRate);
    auto reconstructedNyquist = reconstructSignal(sampledNyquist, nyquistSampleRate, analogSampleRate);
    double nyquistError = calculateError(originalSignal, reconstructedNyquist);

    double subNyquistSampleRate = 1.9 * signalFreq; // 9.5Hz (<2*5Hz)
    auto sampledSubNyquist = sampleSignal(originalSignal, analogSampleRate, subNyquistSampleRate);
    auto reconstructedSubNyquist = reconstructSignal(sampledSubNyquist, subNyquistSampleRate, analogSampleRate);
    double subNyquistError = calculateError(originalSignal, reconstructedSubNyquist);

    std::cout << "Nyquist采样（" << nyquistSampleRate << "Hz）重建误差: " << nyquistError << std::endl;
    std::cout << "亚Nyquist采样（" << subNyquistSampleRate << "Hz）重建误差: " << subNyquistError << std::endl;

    return 0;
}